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ABSTRACT 

A five— parameter gamma distribution (BGD) having two 
shape pareuneters, two location parameters, and a correlation 
parameter is investigated. This general BGD is expressed as 
a double series and as a single series of the modified Bessel 
function. This general BGD reduces to the known special case 
for equal shape parameters. Practical functions for computer 
evaluations for the general BGD and for special cases are 
presented. Applications of the general BGD are to be found 
in reliability theory, signal noise, and meteorology. In 
this paper, applications to wind gust modeling for the ascent 
flight of the Space Shuttle are illustrated. 

This report is from a draft manuscript prepared to meet the 
requirements for a presentation at the American Statistical 
Association (ASA) National Annual Meeting, August 16-19, 1982, 
at Cincinnati, Ohio. The paper for the ASA meeting was ac- 
cepted and will be presented. A version of this report will 
be submitted to the ASA for consideration for publication. 


1 . INTRODUCTION 


Gunst and Webster (1973) discussed the need for a 
practical method of dealing with dependent chl'*sguare vari'- 
ates . They stated that most of the known results discussing 
the bivariate chi-square and related bivariate F-distributions 
were too prohibitive for practical use in that they either 
involved mathematical functions such as Laguerre polynomials 
and cf/nvoluted sums or that they were too restrictive in the 
number of parameters. While Gunst and Webster were primarily 
concerned with dependent chi-square variates, it is obvious 
that the same is true for practical use of dependent geunma 
variates. For example, earlier work by authors such as 
Kibble (1941) , Krishnamoorthy and Parthasarathy (1951) , and 
Down ton (1970) were restricted to the bivariate gamma with 
equal shape parameters. Jensen (1970) and Jensen and Howe 
(1968) were able to consider the unequal shape distributions; 
however, their results are computationally restrictive. 
Recently, McAllister, Lee, and Holland (1981) have been able 
to give efficient methods for computing probabilities from 
Jensen's joint bivariate F-distribution; however, this method 
is still difficult to use. 

In this paper the results given by Gunst and Webster 
have been investigated as a computational model for the 
bivariate gamma distribution with unequal shape parameters. 
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Several functions, such as the conditional distributions and 
probabilities in different geometric regions, are also con- 
sidered. Each of these expressions were of interest in con- 
sidering applications to wind gust modeling for the ascent 
flight of the Space Shuttle. 

Section 2 summarizes the results for the five-parameter 
bivariate density function and considers some special proper- 
ties of the density function. Expressions are derived for 
probabilities over different geometrically shaped regions. 
Some properties of the conditional distribution functions are 
also presented. Section 3 presents an application of the 
bivariate gamma probability functions for wind gust modeling. 
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2. SUMMARY OF FUNCTIONS 


2.1 Bivariate Gamma Probability Density Function 


2.1.1 A Five-Parameter Bivariate Gamma Density Function 


In this section, the results given by Gunst and Webster 
(1970) for the bivariate chi-square are extended to the 
bivariate gamma distribution with unequal shape parameters 
Yl < Y 2 given by 


Yi -1 Y2-I 


f(t^,t2,* Yi/Y2/n) = 


X 


I I 

k =0 j =0 




exp - { (tj^+t 2 )/(l-n) } 


(1-n) ^ r(Y^) rCYj-Y^) 
r(Y2-Yi+k) ^2 


TTy^+J+kT” jl k! 


( 2 . 1 ) 


Where = 3j^x, t 2 = ^ 2 ^' known scale parameters 

a .d n = P'^^Y^Ty^ where p is the correlation coefficient 
between the variables x and y. 


Due to the method of construction by Gunst and Webster, 
the above bivariate gamma density function has five parameters; 
namely, 82 # Yj» Y 2 » n. However, the model is not 

completely general in that the correlation between the 
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variables x and y, denoted by p is restricted to 



Thus, the model is very restrictive whenever Y 2 and p 

is large. To the authors’ knowledge, Moran (1969) has an 
expression for the general five-parameter gamma distribution. 
However, his expression for the density function is compu- 
tationally unfeasible for most practical applications. 


The above function (2.1) can be expressed as a single 
series in terms of the modified Bessel function of the first 
kind, I^(z) as; 


^ y r(Y3-Y^>k) 


Y2-I 

exp -{ (tj^+t 2 )/(l-n) } 


Y 1 

(l-n) T(y^) rCYj'Y^) 

'''’2+k-l ( 



(2.3) 


2.1.2 The Fou ."“Parameter Bivariate Gamma Density Function 

For equal shape parameters, Yi=Y 2 =Y»n=P» and 
2 . 1 becomes : 
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exp -{ (tj^+t2)/(l-p) } 

f(t,,t,; Y,p) » r 

^ ^ (1-p)^ r(Y) 


j=0 


(1-p)^^ r(Y+j) j! 


( 2 . 4 ) 


This function can also be expressed in terms of the Bessel 
function, as: 


f(t-i,t^; Y»P) ~ 


(ti*t 2 ) exp -{ (tj^+t 2 )/(l-p) } 

p(Y"l)/2 (1-p) r(Y) 


X 


^Y-1 \ a-p) J 


(2.5) 


This is the form of the four-parameter bivariate gamma 
density function as derived by Kibble (1941) and reported by 
Downton (1970) . Although no claim can be made that equation 
2.1 is a completely generalized five-parameter bivariate gamma 
density function, it is satisfying to see that 2.1 reduces 
to the known special case. 

Because the four-parameter bivariate gamma density func- 
tion (2.4) is symmetrical with respect to the line tj^ = tj, 
a transformation of the coordinates can be made to take 
advantage of this symmetry as was done by Adelfang and Smith 
(1981) . 
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The expressions to give the finite number of terns suf- 
ficiently large to ensure convergence to an assigned accuracy 
of these infinite series expressions for computational evalu- 
ations can be established. 

2.2 Bivariate Probability Distribution Functions (PDF) 

2.2.1 For the five-parameter density function (2.1) the PDF 
is defined (see Figure 1) as: 



FIG. 1 FSOIAIILITY IN A RECTANGLE 
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After performing the indiceted Integration gives t 


F[ti < tj, tj < tj] “ Pr{tj^ < tj, tj < tj} 


U-n)" V V "^*’'’■'^2-^*''' 

rCY.) r(Y,-YJ ^ ^ jl k! r(Y2+3+ki 

^ ^ k*0 j=0 


) y(Y2+j+J«» r^) • 


Where y(a,x) is the incomplete gamma function which is 


defined as: 


y(a,x) = 


I e- 


dt 


This expression can be expanded in several ways; namely. 


... . V (-1)'' x*‘"” 

= Z ni (a^-nV ' 


y(a,x) = a ^ r(l+a) e ^ 


r(l+a+nl 


P(a,x) = y(a,x)/r(a) . 
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The expression 2.9 is of partioi lar usefulness in ob- 
taining the required integration for sector probabilities of 
Section 2.3. 

2.2.2 For the four-parameter function (2.4), the PDF is: 


PrCtj^ < tj, tj < t*} 


(1-p)^ 

r(Y) 


00 



j=0 


t*/(l-p))y(Y+j^ t*/(l-p)) 

rcY+j) 


( 2 . 11 ) 


Computer programs were developed to evaluate the PDF's 
using the series expressions and by double numerical inte- 
gration of the density functions expressed in terms of the 
Bessel function. Although both methods gave high precision 
results, the evaluations of the PDF expressed as series were 
far more efficient in computational time. 

2.3 Sector Probabilities 

Two expressions are given to determine the probability 
contained in a sector for both the five-parameter function 
(2.1) and the four-parameter function (2.4). 
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The series expression for the integration of Equation 2.12 
is : 


(a»«) 


^2 "'^2 


Y 1 +Y 2 


(l+a) I'<Y^)r(Y 2 -Yi) 


OO 00 00 


y Y V Y r(Y2-Yi+k) r (Y^+Y2+2j+k+n) 

^ ^ ^ j Ik ! \ 2-i+lf +r» »(2. 

r (l+Y2+3+k+n) 


k=0 j=0 n=0 


13) 


where a > 0 o = tan"^ (t*/t*) or a = tan0. 


For the independent case n = 0 in 2.12, the sector 
probability for 2.13 becomes: 


**(g»*) 


(1+g) 


g 

Yl+Yj 


r(Yj) 


^ A (Yj^+Y2'*'*’) 

2 ^ I* (l+Y^+n) 
n=0 ^ 


(l+a)*' 


(2.14) 
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A special case of interest for c«nputer program checkout Is 
to set " 1 in 2.14 and obtain 

^2 

P(a,«) “ [a/(l+ci)l . (2.15) 

In addition, by letting a » tand in 2.14 gives 

P(0,®) = [1 + cotOl ^ (2.16) 

Further, for = 2, 3, 4 ..., the series in 2.14 can also be 
obtained in closed form. 


2.3.1(b) For large a expression 2.12 converges slowly. By 
changing the order of integration of 2.1.2, the probability 
contained in a sector (condition b in Figure 2) is defined 
by: 


P(oi ' ,”) 



f (t^,t2?Yj^»Y2»h)dtj^dt2 


(2.17) 


and the sector probability is obtained as 


(a ' ,™) 


^1 ^2 

g’ ^(l-n) ^ 


d+a’) 


Yi>Y2 


r(Yi)r(Y 2 “ 7 i) 


y Y Y p3‘^*^r(Y2-Yj+k)r(Y^-*‘j)r(Y^+Y2‘*‘2j>k-»‘n)a’^''“" 
k-0 j=o n-0 r ( Y 2 ^ j+k) r ( l+> j+n) ( l^-a ' ) 


(2.18) 
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where a' 



For the independent case, t \ • 0 , equation 2.18 becomes 


P (a ' ,«■) 



d+a’) 


Y1+Y2 


r(Y2> 


00 


I 


r(Yi+Y2+n) a'” 
an.')" 


(2.19) 


Note for the special cases, n = 0; comparing 2.19 with 2.13 
it is seen that the subscripts to the shape parameters are 
interchanged and that a is replaced by a ' . 


2.3.2 Sector Probability for the Four-Parameter Function 


2.3.2(a) For equal shape parameters, Yj^ = Y 2 * Y (see 
Figure 2, condition a), the probability in a sector is; 


P(a 


(X 


w 


( 1 + 0 ) 


(1-p)^ V 

TT7T L 

j=0 


V r(2Y+2j+n) 
Z ]! r (l+Y+3+n) 
n=0 


j+n 


(2.20) 
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and for p - 0 Equation 2.20 bacoaaat 


p . y .r(2Y-»-n).0L! • (2.21) 

(g»-) (1-0)^ rd+Y+nMl+o)” 


2.3.2(b) For equal shape parameters, = Y 2 “ the sector 
probability for condition b shown in Figure 2 is: 


P(a' ,~) 


g’^(l-p)^ 

d+a')^^ r(Y) 


00 00 


x 2 2 


r ( 2 Y^^ 2 i-^•n) 

j! rd+Y+j+n) d+a')^^'*’*' 


( 2 . 22 ) 


2.4 Probability in Triangles 

This section presents expressions to obtain the proba- 
bility contained in a right triangle for two conditions, a 
and b, and for an equilateral right triangle. 
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2.4.1 The Probability in a Right Triangle for the Five- 
Parameter Function 

2.4.1.1(a) The probability contained in a right triangle 
shown as condition a in Figure 3 is obtained by taking the 
definite integral for the second integration of Equation 2.12. 



MS. 3 MSIAIILITY IN A RI8NT TRIAN6LI 

<•) fr (S) <aij* i, <y| 
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2.4.1.1(b) The probability contained in a right triangle 
shown as condition b in Figure 3 is obtained by taking the 
definite integral for the second integration of Equation 2.17. 
The result is; 
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(1-n) 


RT 


d+a') 


Yi+Y 


1-12 


r(Yi) r(Y2“Yi) 


00 OO 00 


X I I I 

)c=0 j»0 n=0 


r(Yi-»-j) 


, j+lc+n 


jl k! r(Y2+j+l^) r(l+Yi+j+n) 


X y(Yi+Y2+2j+k+n, (^) t*) 


(2.25) 


For ri = 0 Equation 2.25 becomes: 


a 


Jl 


RT 


( I+a ' ) 


Yi+Y 


1^*2 


^ yCYi+Yj'*'”' 

Z rd+Yi+n) 

n=0 


a 




(2.26) 


d+a') 


n 


2. 4. 1.2 The Probability in a Right Triangle for the Four- 
Parameter Function 


2.4.1.2(a) The probability contained in a right triangle 
for the equal shape parameter bivariate gamma function for 
condition a shown in Figure 3 is; 


For Yi = Y, = Y Equation 2.23 reduces to; 
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p 


RT 
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(1+a)^^ r(Y) 


V V p’ 

r(i+Y+j+») 

j-0 n»0 



X y(2Y+2j+n, (^) tj ) 


( 2 . 27 ) 


and for p * 0 Equation 2.27 becomes; 


P 


RT 


a 




(1+a) 


(Y) 


V 

n*0 


y(2Y+n, (1+a) t*)a” 
r(l+Y) (l>a)" 


( 2 . 2 ») 


4.4.1.2(b) For condition b shown in Figure 3, the probability 
in a right triangle is: 


For Yj^ = Y 2 “ 7 


oo Ob 


»•’' (i-p) “ "V ^ 


„.j*n 


y(2,+21.n, [^) tj) 


(2.29) 
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emd for p - 0 Equation 2.29 becomes 


\ Y{2y+n, (1+a') tj ) 

(l+a’)^^r(Y) m+Y+nl 


X — ' 


(1+a') 


2.4.2 For Equilateral Right Triangle 
2. 4. 2.1 For the Five-Parameter Function 


(2.30) 


F The probability contained in an equilateral right tri 

I angle shown in Figure 4 for the five-parameter function is 

obtained from (2.31). 



Wm' . 







ERT 


OWOINM. ® 
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I I fCtj^rtj; YirY2»n)dt2dtj^ 


The indicated integrations of 2,31 yield: 


ERT 


Y *^Y 

t? ^ ^ exp -{t*/(l-n)} 


( 1 -n) ^ r(Y 2 -Yi) r(Yj^) 



Y Y Y r(Y 2~Y]^+h) r(Yj^-fj) r (i+Y2^.j-*-ig-»-n) 

J.t'o n =0 jl hi r(l+Yi+Y 2 + 2 j+l«+n) 


32 ) 


For n = 0 Equation 2.32 reduces to: 


Yt+Y 2 -t? 

^ERT= ^1 


£ Y r(i+Y2+>^)ti 
2 ^ TTI+Y^TY^+ny 


( 2 . 33 ) 


n =0 


2 . 4 . 2. 2 For the Four- Parameter Function 

For equal shape parameters, Equation 2.32 reduces to: 
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For = Y2 “ 


exp ~{t£/(l-p)} 
(1-p)^ r(Y) 


ERT 


00 00 


-I S 

j=0 n=0 


r(Y + j) r(i+Y+ j+n) tj 


2 j+n 


jl r(l+2Y+2j+n) 


(2.34) 


a.nd for p ^ C ^ Ecjustion 2 * 34 b©con\©s : 


2 ^ r(l+Y+n)t*’^ 

^ERT = Z r(l+2Y+^r 

n=0 


(2.35) 


2.5 Properties of the Conditional Distribution 
2.5.1 Conditional Probabilities 

Because the five-parameter 'jairana function 2.1 is not a 
sytranetrical function, there are two expressions for the con- 
ditional distribution. They are for (tj^ 1 tj) and (t 2 1 tj^) 
given by 2.36 and 2.38. 
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r(Y 2 ^ exp -\ t =^\ 

r(Yi) rcYj-Yj^) 



j 


y(Yi+jr Vd-n) 

J 1 



r (Y 2 "Yi+k) 
r(Y 2 +j+Jt) ki 


(2.36) 


By taking the limit of the given value in 2.36 as t -2 0 

produces : 


Pr{tj^ 1 I ^2 ^ 


y(Yi, tj^/d-n)) 

nrp 


(2.37) 


Pr{t 2 < t 2 


~ L J ; Yj^»Y 2 » h) 


d-n) 


Yo-Y 






1 

JT 


^ ^n’' r(Yj-Yj.k) 

Z k! r(Y,+i+k) 

k=0 ^ 


y(Y2+j+k, t2/d-n)) I • 


j 


(2.38) 
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and upon taking the limit as ® 2.38 gives 


Yj-l 


Pr(tj < 1 » 0) - 


^ k r (Yo-Yi+Jt) 

2 ITT r{Y,+k) y(Y 2 +l^» t 2 /(l-n)) . 


k=0 


(2.39) 


For the equal shape parameter bivariate gamma function, 
the conditional PDF is 


Pr{t 2 < t 2 1 tj^ = tj; Y,P) = exp 




X 


2 

K = 0 


(1-p) 


It ^ ^ 

^ y(Y+k, t2/(l-p)) , 


k k! 


(2.40) 


and for the limit as tj -*• 0 


Pr{t2 1 ^2 1 -*• 0; Y,p) » 


Y-1 


(1-p)^ r(Y) 


X exp - { t 2 / (1-p) ) 


(2.41) 


Because this is a symmetrical function, the conditional 
distribution for { tj^ 1 t 2 > can be obtained by interchanging 
the subscripts in Equations 2.40 and 2.41. 
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Equations 2.37, 2.39, and 2.41 are interesting results. 
This is because the original bivariate gaiwia functions (2.1 
and 2.4) are bounded at the origin and yet the conditional 
functions are not bounded at the origin as the limit of the 
given values approaches zero. This property of the condi- 
tional functions must be reasoned with when they are used as 
models for physical data. 

2.5.2 Moments of the Conditional Function 

For the four-parameter (equal shape parameter) bivariate 
gamma density, f(tj^, t 2 i Y, P) , the conditional moments are 
sufficiently simple to be derived. For this case, the non- 
central conditional moments are expressed as 




r(ytj) (1-p)^ 

r(Y) 


-(pt*)/(l-p) 

e 


X j^F^ (Y+j ;y; (ptj)/(l-p) ) , 


(2.42) 


where |Fj^(a ; c ; z) is the confluent hypergeometric function 
of Hummer's function as defined by W. Magnus et al (1966). 

From 2.42 the conditional mean and variance are given 
by 

E(^_^ : - t*) =• 1-Ph + PtJ/(l-»')| (2.43) 
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Var (tj I = t|) - (l-p)2 [y + 2pt{/(l-p)] 


(2.44) 


Note that the conditional moments exist when the given value 
t| approaches zero in Equations 2.43 and 2.44. Also, it is 
observed from 2.43 and 2.44 that the regression of tjOnt^^is 
not homogeneous over t^^. 

2.6 Illustrations and Properties of the Bivariate Gamma 

Distribution 

This section presents graphical displays of the bivariate 
gamma density function to illustrate some of the properties 
of this function. 

Figures 5a, 5b, 6a, 6b, 7a, and 7b are schematic illus- 
trations of the contours of equal probability density to 
show the variations with the shape parameter, y, and the 
correlation coefficient, p, for the four-parameter function, 
f(ti, tj,* Y. P) » Equation 2.4. In all Figures 5a-7b, the 
density contours are symmetrical with respect to the line 
ti = t2* 

For Y < 1, Figures 5a and 5b, the density is a maximum 
for small values of t^^ and and becomes undefined at 

. tj . 0. 

For Y ® 1, figures 6a and 6b, which is the bivariate 
exponential density function, the maximum density is at 
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t * t » 0. For Y * p * 0 (Figure 6a ) $ the contours 
for equal probability density form straight lines normal to 


the line t^^ = t 2 * 


For Y > 1 (Figures 7a and 7b) , the mode exists and 
lies on the line t^^ = t 2 * (Additional characteristics of 
the mode and the density of the mode are shown in Figures 10 

and 11.) 

Because the equal shape parameter bivariate gamma density 
function 2.1 is symmetrical with respect to the line = t 2 » 
there is some analytical advantage in expressing this func- 
tion in a coordinate system rotated by 45 degrees. In terms 
of the Bessel fu. ct ion. Equation 2.5 in the rotated coordinate 

system is 

Yzi 

(z ^-z ^)^ 

‘ (l-p)(2p)''-^'^' r(y) 


X exp 


- {/2 z,/(l-p)} 


TT^p) 


(2.45) 


where ~ ~T ^®x ^ 


,2 = ^ (By Y - 6^ X) 


and > 0, Z 2 1 0, z^ > Z 2 ; v > 0 and 0 < p < 1. 

This Bessel function is available as a computer subrou* 
tine to aid in evaluating 2.45. 
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Computer graphics for the evaluation of 2.45 for con- 
tours of equal probability density are shown in Figures 8a » 
8b, and 8c for y ® 2 with p* 0.25, 0.50, and 0.75; and in 
Figures 9a, 9b, and 9c for y = 3 with p * 0.25, 0.50, and 
0.75. The density contours were obtained by the following 
steps : 


1 . 


2 . 


3. 


Iterative evaluations of 2.45 with S 2 * 0 to 
obtain the model value, 

Dividing the into 11 increments, Az^^, and 
evaluation of 2.45 with Z 2 = 0 and Az^^ to obtain 
the 11 density values to be contoured. 

The contouring for the 11 density values was by 
iterative evaluations of 2.45 for the valid range 
of and Z 2 for each contour. 


Interesting special cases of 2.45 can be expressed. One 
is for Y = 1.5. For y = 1.5, the Bessel function reduces to 
a hyperbolic function and the PDF becomes: 


f(z,,z-; Y = 1.5, p) = 


nip(l-p) 1 


X exp 




sinh 


1 


/7B (z^2^Z2^)** 
— 


(2.46) 
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from which the inoder z-^, is 


2 = f tanh ^ {/p} 

1 /Ip 


and the modal density, D, is obtained in closed form as; 


D = 


2 exp -{xT^ ^l} 


n/p (1-p) 


{li^ ^i} 


The values of z~ ^ function of z^, for p fixed and density, 


D = constant, is expressed as: 


, 2 _ , 2 q-p) 
^2 "1 


sinh ^ /p (1-p) • D • exp|^_p'y z-^ 


By using the identity for the inverse hyperbolic sine 


function, 2.49 becomes 


2 _ , 2 (1-p) 

*2 " ^1 JF" 


In \C + + 1 


where C = ^ /pTT-p) • D • exp z^} 


Then Z 2 = + * ^2 
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Equation 2.49 or 2.50 can be used as a special case to cosi- 
pute contours of equal probability density. 

Another special case is to let p * 0 and for y > 1 in 
2.45 to get 




1. P 


0) 




- *2^ 


Y-1 


2^"^ [r(Y)] 


exp -{/7 
5 


(2.53) 


and 

f{z^,Z2 = 0, Y > 1, P 


0) 


exp - {/5 
2^"^ ir(Y)]^ 


from which the mode is 
= /2 (Y-1) 


(2.54) 


(2.55) 


and density at the mode is 


/ tv2(y-1) 

5 = X — exp - {2 (y“D) 

(r(Y)]'^ 


(2.56) 


Then 


2 ,2 

Z2 = 


as a function of Zj^ and fixed density, D, is 


2^"^ (r(Y)l^ * D • exp {/2 z^} 


1 
Y- 


(2.57) 
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The location of the mode and the density value at the mode 
are informative properties of skewed and usually bounded 
probability distributions. 

For the four-parameter gamma probability function (2.45), 
the location of the mode, versus y for correlations, p, is 
shown in Figure 10. Here it is observed that the mode, z^, 
increases with increasing y and decreases for increasing p. 

The value of the density at the mode (Figure 11) decreases 
with increasing y and increases with increasing p. 

Three-dimensional computer graphics are shown in 
Appendix A for equal and unequal shape parameter gamma densi- 
ties as the correlation varies . Also illustrated are contours 
for equal probability densities. Although the abscissa and 
ordinate scales are unequal, it can be seen that more density 
is above the line, tj^=t 2 , than below this line for Y 2 > y^. 

Within Section 2 there are elements on the properties 
of the bivariate gamma probability distribution which are 
believed to be original contributions to the better under- 
standing of this function. Just as it is important to under- 
stand the characteristics of a data sample to make statistical 
infevcncos and to model the sample by a probability function, 
it is also impt)rtant to understand the properties of the 
probability function in the selection of other alternatives 
as a model for the Scunple data. 
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3 . APPLICATION 


Aerospace launch vehicles experience structural loads 
that are imposed by wind variability along the vehicle tra- 
jectory. The largest loads are often associated with a 
discrete gust that stands out above the general disturbance 
level (Jones, 1971). To ensure a margin of safety, informa- 
tion about the amplitude and length scale of discrete gusts 
and how they vary with altitude and season is required. A 
previous discrete gust model used in aerospace engineering 
applications proposes a gust with a quasi-square wave shape, 
a variable length scale, and constant amplitude (Daniels, 

1973). Other studies (Adelfang, 1970; Fichtl and Perlmutter, 
1976) based on analysis of detailed Jimsphere wind profiles 
between the surface and 16 km have concluded that gust ampli- 
tude increases with altitude above 9 km. A new model described 
herein treats gust amplitude and length scale as the variables 
of a bivariate gamma distribution. Parameters of the model at 
various altitudes are estimated from samples of gust data de- 
rived from wind profiles. 

3.1 Data 

Detailed wind profiles at Cape Canaveral, Florida, mea- 
sured from the earth's surface to 18 km at 25-meter altitude 
intervals with the Jimsphere FPS-16 measurement system (Camp, 
1971) were used. The sample consists of 150 profiles per 
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month for the months of February, April, and July. Data 
suitable for analysis were derived from these profiles by 
application of three band-pass Martin-Graham digital filters 
(DeMandel and Krivo, 1971) . The filtered profiles defined 
here as residual profiles contain wavelengths within the 
90-420, 420-2470, and 2470-6000 meter bands. A set of zonal 
wind component, u', residual profiles calculated from a 
Jimsphere profile at Cape Canaveral (2/23/71, 1445 GMT) is 
illustrated in Figure 12. 

The definition of gust used in this study satisfies the 
objective to provide data suitable for detailed analysis of 
singularities or discrete perturbations that are often ob- 
served in Jimsphere wind profiles. According to the conven- 
tional approach, all the amplitudes in the filtered profile 
would be defined as gusts. In this study, only the largest 
residual with equivalent sign to the residual at reference 
height, Hq, is defined as a gust; as illustrated in Figure 13, 
the gust, u', occurs between successive zero crossings at 
and H 2 ; L^, is defined as the gust length. 

3.2 Results 

3.2.1 Marginal and Conditional Distributions 

Examples of expected (gamma) and observed PDF’s for absolute 
gust component amplitude, lu'|, and associated gust length. 
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Luif are illustrated in Figures 14 and 15; the garnna PDF's, 
illustrated with solid lines, were calculated from aMonpla 
_ estimates of the parameters y and 6 listed in the figure 

^ insets. The relatively small gust amplitudes in the 90-420 

- meter wavelength band and the good agreement between gamma 

^ and observed PDF's are clearly indicated. 





I 
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A further indication of goodness of fit was obtained 
from calculations of the X statistic from expected and ob- 
served populations. Although the outcome can be affected by 
the arbitrary choice of the size and number of class intervals, 
the results to date have supported the conclusion that the 
marginal distributions are gamma distributed. A partial 
summary of the results is illustrated in Table 1. Results 
for other months and other wavelength bands have yielded a 
similar conclusion. 


I 


The variation of gust component amplitude with altitude 
and season is an important consideration in aerospace design 
and operations. Illustrations of how the univariate gamma 
distribution is used to describe this variability are given 
in Figures 16 and 17. During February at Cape Canaveral, 
u component gust amplitude increases with altitude between 
4 and 12 km; the variation is relatively small between 4 and 
8 km and large between 10 and 12 km (Figure 16) . As indi- 
cated in Figure 17 , gust amplitude is much smaller in summer 
(July) and early fall (October) than in Winter (February). 
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L^,, are illustrated in Figures 14 and 15; the gamma PDF's, 

illustrated with solid lines, were calculated from sample 

estimates of the parameters y and $ listed in the figure 

insets. The relatively small gust amplitudes in the 90-420 ! 

meter wavelength band and the good agreement between gamma 

and observed PDF's are clearly indicated. 

A further indication of goodness of fit was obtained 

2 

from calculations of the X statistic from expected and ob- 
served populations. Although the outcome can be affected by 
the arbitrary choice of the size and number of class intervals, 
the results to date have supported the conclusion that the 
marginal distributions are gamma distributed. A partial 
summary of the results is illustrated in Table 1. Results 
for other months and other wavelength bands have yielded a 
similar conclusion. 

The variation of gust component amplitude with altitude 
and season is an important consideration in aerospace design 
and operations. Illustrations of how the univariate gamma 
distribution is used to describe this variability are given 
in Figures 16 and 17. During February at Cape Canaveral, 
u component gust amplitude increases with altitude between 
4 and 12 km; the variation is relatively small between 4 and 
8 km and large between 10 and 12 km (Figure 16) . As indi- 
cated in Figure 17, gust amplitude is much smaller in summer 
(July) and early fall (October) than in Winter (February) . 
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FEB, 12 KM, LU 

PARAMETERS OF THEORETICAL 
DISTRIBUTION 

WAVELENGTH ? f 

BAND 

(hi) (1/m) 

90-420 2.781 .026B 

420 - 2470 4.199 .00601 

2470 - 6000 4.757 .00151 
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FIGURE 15. GAMMA AND OBSERVED PDF OF GUST LENGTH, 

Lu'. FROM 90-420m, 420-2470m, 2470-6000m WAVaENGTH BAND 
FILTERS AT 12km ALTITUDE, FEBRUARY. CAPE CANAVERAL, aORIDA 





ALTITUDE 12km 




figure 16. GAMMA PDF OF GUST COMPONENT AMPLITUDE |u‘| , FROM 420’2470 m 
WAVELENGTH BAND FILTER AT 4, 6 ‘-14km ALTITUDES, FEBRUARY, CAPE CANAVERAL, aORlOA 
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FIGURE 17. PDF OF GUST COMPONENT AMPLITUDE, ju'l, FROM 420-2470m 
WAVELENGTH BAND AT 12km ALTITUDE FOR FEBRUARY, JULY, OCTOBER, CAPE CANAVERAL, 
FLORIDA. 








Tabid 1. Rdsulbs of Tdsbin^ the Hypothesis that the Samples^ 
of u and v Component Absolute Gust and Gust Length 
are Drawn from Gamma Distributed Populations 


Wavelength 

Band 

(m) 

A/R* 

lu'l 

|V| 


^v' 

All 

Variables 

90-420 

4/2 

6/0 

5/1 

5/1 

20/4 

420-2470 

5/1 

6/0 

6/0 

6/0 

23/1 

2470-6000^ 

4/1 

5/0 

4/1 

4/1 

17/3 

All Filters 

13/4 

17/0 

15/2 

15/2 

40/8 


‘Samples for month of February at six altitudes (4, 6, 

14 km) for Cape Canaveral, Florida. 

* Ratio of cases accepted to cases rejected utilizing the 
test at the .05 level of significance. 

3 

Five altitudes (6, 8, ..., 14 km) 



Aerospace vehicle control systems can respond critically 


to large wind perturbations at certain wavelengths. The con- 
ditional PDF of gust component amplitude, |u'|, given a gust 
length, for various altitudes and months is useful in 

vehicle design. Conditional PDF's calculated from Equation 
2.36 using sample estimates of the parameters from Cape Canaveral 


band pass (420-2470) gust data at 12 km are illustrated in 
Figure 18. Percentiles of the conditional PDF plotted as a 
function of the given gust length are illustrated in Figure 19. 
It is indicated that each percentile can be approximated by an 


®mpitical linear function of for this example. 
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CAPE CANAVERAL, aORIDA. 




99th PERCENTILE 


OWGINAL ^ 

OF POOR QUALITY 



figure w. percentiles of conditional gust component amplitude, |u'1 , GIVEN GUST 

LENGTH. L.j„ FROM 420 - 2470m WAVELENGTH BAND FILTER AT 12km ALTITUDE FEBRUARY. 

CAPE CANAVERAL, FLORIDA 


3.2.2 Bivariate Distributions 

An observed bivariate distribution of gust varied>les is 
illustrated in Figure 20. The nondimensional variables 
and ?2 are calculated from u component absolute gust, |u'l, 
and associated gust length, L^,, data at 12 km over Cape 
Canaveral for a wavelength range of 420-2470 meters. The 
probability contained within rectangular domains was calculated 
according to Equation 2.6. Rectangular domains are well suited 
for preliminary evaluation of the goodness of fit of the model 
to the observed distribution. A comparison is given in 
Table 2 of the observed number of occurrences (counted within 
2x2 cells in Figure 20) and the expected number if T^^ and T 2 
are bivariate gamma distributed. There is good eigxeement 
between observed and expected cell counts. 



The domain of bivariate gamma distributed variables can 
be divided into a number of equal-area sectors. The probability 
contained within each sector is obtained by calculating the 
difference in probability between two partially overlapping 
sector probabilities. For sectors of equal area, the vari- 
ation of probability from sector to sector is a measure of the 
relationship between the shape parameters of the marginal 
distributions. Equivalence of shape parameters yields a sym- 
metry in the variation of the sector probabilities! whereas, 
nonequivalence yields asymmetry. The change in symmetry is 
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'able 2. Observed* and Expected* Number of Occurrences within 2x2 Cells for 



Cape Canaveral at 12 km during February for a wavelength range of 420-2470 meters; 
integer values in upper left of table entries for each cell. 

Expected for correlated bivariate gamma distributed variables; values in lower 
right of table entries (calculated from Equation 2.6). 

T, = B,|u'I, T, = 6-> L refer to Figure 20 caption for parameter estimates. 




s 

clearly indicated in Table 3. A case illustrating asymmetry 
in an observed distribution is illustrated in Figure 20; the 
significant difference in the sample estimates of the y param- 
eters (listed in the figure inset) is consistent with the 

; asymmetry of the plotted data. Nearly 90 percent of the data 

< « 

are within the lower half of the domain illustrated in 
Figure 20. 


Table 3. Sector Probabilities for Equal Area Sectors 
(Calculated from Equation 2.14), Yj^ = 2 


Sector Range 
Deg. 

2 

^2 

3 

4 

5 

0-15 

.11510 

.03177 

.00829 

.00208 

15-30 

.18875 

.11054 

.05518 

.02531 

30-45 

.19615 

.17020 

.12403 

.08198 

45-60 

.19615 

.22211 

.21056 

.18047 

60-75 

.18875 

.26696 

.31588 

.33770 

75-90 

.11510 

.19843 

.28606 

.37245 
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4 . CONCLUDING REMARKS 


Although the main topics stated in the abstract and 
introduction have been investigated, the*:e remains a need 
for further investigations to address such topics as param- 
eter estimation and hypothesis testing for data sas^lea taken 
from correlated bivariate gamma distributed variables. 

Three commonly used procedures for estimating the 
parameters for the univariate gamma distribution are: 

1. The moment method, 

2. A maximum likelihood method by Thom (1966), and 

3. A polynomial approximation given by Greenwood and 
Durand (1960) . 

Bury (1975) gives a detailed discussion on parameter estima- 
tion in the gamma distribution. 

Very little can be found in the open literature concern- 
ing hypothesis testing for the parameters of the gamma distri- 
bution. In roost cases, the results given in the literature 
assume that the shape parameters are known and tested for the 
scale parameter. In this report, a test for equality of 
shape parameters in the presence of correlation was needed. 
However, this problem has not been satisfactorily resolved 
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^nd remains unresolved to the best of our knowledge. That 
is, there is not a test for 




whenever (x,y) % T p) and p ft 0. Moran (1969) 

A X y y 

has given an expression for the joint likelihood function and 
has indicated a procedure for testing Hq ' ~ ' 

Further discussions on the properties of the functions 
presented in Section 2 and additional examples in Section 3 
may be helpful to those investigators who are interested in 
making applications of these distribution functions. 
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APPENDIX 


GRAPHIC ILLUSTRATIONS 


This appendix presents three-dimensional computer 
graphics of bivariate gamma density (Figures A1-A17) and con 
tours of equal probability density (Figures A18-A34) for 
equal shape parauneters (2.4) and unequal shape parameters 
(2.1) as the correlation varies. 
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figure A1& BIVARIATE GAMMA DENSITY CONTOURS, 
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FIGURE A19. BIVARIATE GAMMA DENSITY CONTOURS, 7. - % ' 1- 5, p- a 25 
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FIGURE A2l BIVARIATE GAMMA DENSITY CONTOURS. 7i - T 2 • 1 5. p- a 75 
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figure A25l bivariate gamma density contours, 7j- 72-3.0, p-a75 
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FIGURE A27. BIVARIATE GAMMA DENSITY CONTOURS, Ti'lO, To* 0.0 
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FIGURE A3a BIVARIATE GAMMA DENSITY CONTOURS, Tj ■ 1, 7, ’ 3, f}’ ft 75 
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FIGURE A34. BIVARIATE GAMMA DENSITY CONTOURS 
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